Analysis of wasp-waisted hysteresis loops in magnetic rocks 
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The random-field Ising model of hysteresis is generalized to dilute magnets and solved on a Bethe 
lattice. Exact expressions for the major and minor hysteresis loops are obtained. In the strongly 
dilute limit the model provides a simple and useful understanding of the shapes of hysteresis loops 
in magnetic rock samples. 
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, Hysteresis is a non-equilibrium effect |f| . If a system is driven by a cyclic force that changes faster than the system 
■ can adjust to it, then the response does not move up and down on a single path but rather makes a hysteresis loop. 
Theoretically the area of the hysteresis loop should vanish as the frequency of the driving field goes to zero but several 
systems with quenched disorder show large hysteresis even at the slowest driving rate. This behavior arises from the 
presence of numerous mctastable states in the system that are separated from each other by energy barriers much 
T , larger than the thermal energy The metastablc states correspond to local minima in the free-energy landscape of the 
system. The system remains practically trapped in a local minimum and is unable to attain thermal equilibrium over 
observation times. However, it can be made to jump from one local minimum to another if a sufficiently strong force is 
applied to it. We focus on magnetic systems. The magnetization induced by a cyclic field traces a hysteresis loop. The 
loop is essentially a locus of magnetizations of mctastable states along the trajectory. Its shape and area are clearly 
objects of practical interest because these determine the rate of energy dissipation in the system. Somewhat less 
1 , ' obvious but equally important is the fact that the hysteresis loop also contains information regarding the distribution 
of local free-energy minima and the energy barriers between them. Sethna et al introduced the non-equilibrium 
q random-field Ising model of hysteresis that not only reproduces the shapes of hysteresis loops 'pleasantly familiar to 
, the experimentalist' but also provides an understanding of other aspects associated with it e.g. Barkhauscn noise, 
l— "™ 1 ■ return point memory, and critical-point phenomena 0, Q ■ 

I | Originally the random-field Ising model was introduced to study the effect of quenched positional disorder on the 
critical behavior of a system in thermal equilibrium [3] . It showed that even an arbitrarily small amount of disorder 
I/"") ' raises the lower critical dimension of a system. The lower critical dimension is the dimension below which a system 
^sO , can not possess long-range order in a state of thermal equilibrium. Above its lower critical dimension, it may evolve 
' into an ordered state at low temperatures but its approach to thermal equilibrium is not smooth. This is due to the 
presence of a large number of local minima in the free-energy landscape of disordered systems [f| • The local minima 
are surrounded by high barriers. The barrier heights are random but much higher than the thermal energy of the 
' system. Hence the approach to thermal equilibrium is very slow and sporadic. Indeed the system may not reach 
equilibrium over practical time scales. This makes the determination of the equilibrium state (the global minimum of 
free-energy) a difficult task numerically as well as analytically. The dynamics of the approach to thermal equilibrium 
in the random-field Ising model has been well studied @, but the progress has been slow due to the difficulty of 
the problem. 

The difficulty of thermal equilibration is sidelined if the model is used to study hysteresis. Practically speaking, 
even in the zero frequency limit of the driving field, hysteresis is observed over time periods that are much shorter than 
the time required for the system to equilibrate. Thus the thermal relaxation process and the global minimum of the 
free-energy are not of primary importance in this case. In the non-equilibrium random-field Ising model proposed by 
Sethna et al , the time required by the system to equilibrate is set equal to infinity. Equivalently, the temperature 
of the system is set equal to zero. The metastable states under the stochastic thermal dynamics thus become stable 
states under the athermal zero-temperature deterministic dynamics. This does not compromise with the essential 
physics of the problem. There is an argument based on the renormalization group theory that the phase transition 
in the equilibrium random- field Ising model is controlled by a stable zero-temperature fixed point [7[. Several key 
features of hysteresis observed in disordered fcrromagnets as well other systems whose dynamics is characterized by 
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avalanches are very well reproduced qualitatively and even quantitatively by the zero-temperature non-equilibrium 
random field Ising model [8(. The deterministic dynamics also makes the model amenable to an exact solution in 
some special cases . 

This paper generalizes the non-equilibrium random-field model of hysteresis to dilute magnetic systems. We solve 
the dilute version of the model exactly on a Bethe lattice and apply the results to explain the shapes of hysteresis loops 
of magnetic rocks. We choose magnetic rocks for two reasons: (i) these are natural realizations of very dilute magnetic 
systems, and (ii) have not received the same attention in the physics literature as in geolo gy. It is recognized that 
rock magnetism arises from a few percent or less of magnetic minerals present in the rocks [l2l [l3j . The commonly 
occurring magnetic minerals are magnetite (FesO^), maghemitc (7^203), titanomagnctitc (F e2- y Ti y O^) 1 pyrrhotite 
(FerSs), greigite (Fe 3 Si), hematite (aFe 2 3 ), and goethite (FeOOH). The composition [14j is generally determined 
by breaking the rock sample. Recently, hysteresis measurements have been employed as a a possible nondestructive 
alternative. The motivation for this comes from the fact that the hysteresis loop of a rock sample as a whole is 
generally different from that of its magnetic constituents in their pure form. Pure magnetic minerals are mostly 
ferromagnetic or ferrimagnetic and therefore their hysteresis loops are similar to those of iron (Fe). On the other 
hand, hysteresis loop of a rock sample as a whole can have unusual shapes including a wasp-waisted shape that is 
constricted in the middle The wasp-waisted shape is thought to arise from the fact that the grains dispersed 

in the rock have a distribution of sizes, domain states, and coercive fields [l5l - [27j . One would like to determine the 
distribution of magnetic grains from the hysteresis loop of the rock sample. It is not immediately obvious if this is 
feasible. A given distribution of grains and a set of interactions between them will produce a unique hysteresis loop. 
However, the inverse relationship of the hysteresis loop to the contents of the rock is not necessarily unique. This 
is because the detailed information concerning the grains has to be integrated over in order to obtain the hysteresis 
loop. Nevertheless there are practical advantages in exploring the extent to which we can deduce the composition of 
a rock from its hysteresis loop. As a first step towards this goal, one would like to study simple models of hysteresis 
and explore how their predictions vary with the parameters of the model. 

There are a number of models of hysteresis in magnetic rocks in the field of geology. The majority of these use 
qualitatively similar assumptions and differ from each other only in detail. As an illustration we consider the model 



studied in reference [17( . Grains of magnetic minerals of various sizes are assumed to be randomly dispersed in the 



rock. The grains are frozen in their positions but their magnetization can rotate or flip under a driving field as well 
as under thermal fluctuations. Interactions between the grains are neglected. This is presumably justified because 
the grains occupy only a few percent of the total volume of the sample and therefore they may be well separated from 
each other to interact significantly. The size of a grain determines the quality of its response to the applied field. 
Small grains (say 30 nm or less) behave like a paramagnet and respond to a cyclic field without hysteresis. Larger 
grains behave like a ferromagnet and respond to the cyclic field with hysteresis. This is because the size of the grain 
is related to the size of the energy barriers that stand in the way of smooth rotation or flipping of magnetization 
within a grain. These barriers are small for small grains and large for large grains. The thermal energy of the system 
provides the criterion for deciding whether a grain is small or large. If the barriers are small in comparison with the 
thermal energy, the magnetization is able to rotate freely under thermal fluctuations and is able to attain thermal 
equilibrium. In this situation the average magnetization is zero in zero applied field. In a non-zero applied field it is 
given by the Langevin function for a paramagnet. Larger grains are not able to attain thermal equilibrium over time 
scales of the experiment. Their response to the cyclic field is the non-equilibrium response in the form of a hysteresis 
loop. A hysteresis loop may be characterized by a number of parameters such as the saturation magnetization, the 
remanent magnetization, and the coercive field. These parameters are related to the size, shape and material of the 
grain. Some variants of the model of hysteresis in a rock sample consider directly a distribution of the coercivities 
of the grains, others a distribution of the sizes that is translated into coercivities through an assumed relationship 
between the two. The task of the models is to reproduce the experimentally observed shapes of hysteresis loops for 
a reasonable choice of the parameters of the model. Some models achieve this task by considering grains of different 
sizes [13, HH, others b y s trongly contrasting coercivities [ill [llj]. A few models have also included interactions 
between the grains [23l425| . A number of models reproduce major hysteresis loops similar to those observed in the 
laboratory experiments. Thus the difficulty is not that we do not have a model of rock magnetism but that we have 
several. The question naturally arises if we can distinguish between these models. Some authors have suggested that 
the comparison of experimental first-order reversal curves (minor hysteresis loops) with the predictions of different 
models may serve to distinguish between them (24[. 

As stated earlier, we adapt the non-equilibrium random-field model of hysteresis to dilute magnets. It offers 
a framework for understanding hysteresis loops of magnetic rocks and the relationship between different models 
employed for the purpose in the field of geology. The non-equilibrium random field Ising model 0, Q is defined on a 
lattice whose sites are occupied by an Ising spin (a binary variable that represents a unit domain with magnetization 
±1). The magnetization of a unit domain is allowed to flip (up/down) rather than rotate continuously. Each spin 
interacts with its nearest neighbors and experiences a uniform external field as well as a Gaussian quenched random- 
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field with average value zero and standard deviation a. This model along with the zero-temperature Glauber dynamics 
has played a key role in understanding several aspects of ferromagnetic hysteresis including Barkhausen noise, return 
point memory, and scale invariant avalanches characterizing critical hysteresis. The ingredient we add to this model 
is the random dilution of magnetic sites. We are interested in the limit of large dilution when only a few percent of 
the sites are occupied by spins. In this case the spins form small isolated clusters of different sizes. The absence of a 
large spanning cluster on the lattice precludes scale invariant avalanches or critical hysteresis but reproduces shapes 
of hysteresis loops that are commonly seen in magnetic rock samples. The size distribution of scattered clusters on 
the lattice- is determined by the random occupancy of the lattice sites by spins. A cluster may be thought of as a 
magnetic grain in other models of rock magnetism but here there is no need to make a separate assumption for the 
distribution of grain size. It is also easily understood why small clusters behave somewhat like paramagnets and 
larger ones like ferromagnets. Spins flip up or down whenever the net field at their site changes sign. An isolated spin 
(smallest cluster) has no memory and behaves like a perfect paramagnet. A spin connected to other spins is influenced 
by them in addition to the on-site magnetic field and therefore it shows hysteresis. The random field Ising model of 
ferromagnetic hysteresis can be solved exactly on a Bethe lattice and gives important insights into the behavior of 
the model 0, [13]. In the following we extend this solution to the dilute case. The exact solution in the dilute case is 
convenient in studying the effect of changing various parameters on the shape of hysteresis loops without performing 
time consuming simulations of the model. 



II. THE MODEL 



The dilute random-field Ising model is defined by the Hamiltonian, 



H = -j'^CtCjSiSj - hjCjSi - h'^c i S ll (1) 

i,j i i 

where J ( J > 0) is a ferromagnetic exchange interaction and the sum is over nearest neighbor sites {i,j} of a lattice. 
The restriction to nearest neighbor interactions means that the model applies to systems in which long range dipolar 
interactions are negligible. Si = ±1 is an Ising spin, hi is a random- field, and h is a uniform external field. The 
random-field hi is drawn from a Gaussian distribution with mean value zero and standard deviation a. The quantity 
Ci is a random variable taking the value 1 with probability c and with probability 1 — c. Thus c is the concentration 
of lattice sites occupied by spins. The quantities {hi} and {c{\ are quenched and therefore remain unchanged under 
the evolution of the system. The spins are the dynamical variables. These are governed by the zero-temperature 
Glauber dynamics at discrete time steps t, 



S t (t + l)= sign 



J CjSj(t) + hi 



(2) 



The sum on the right-hand-side is over nearest neighbors of site i. At a fixed applied field h, this dynamics lowers 
the energy of the system iteratively and takes it to a stable fixed-point {S*(h)} such that for each lattice site i 



S*( h ) = si g n 



jJ2 c ^;(h) 



(3) 



We are interested in the hysteresis loop when the applied field h is cycled from — oo to oo and back to — oo 
infinitely slowly. Numerically, we start with a sufficiently large and negative h, so that the system is at the fixed-point 
{S*(h = — oo) = —1}. Now h is increased slowly until this fixed point becomes unstable. We hold h fixed at this value 
and allow the system to evolve under the iterative dynamics until it reaches a new fixed point. The new fixed-point 
is characterized by its magnetization per site m*(h), 



m*(h) = -^CiS:(h) (4) 

i 

The process is repeated i.e. we start with a fixed-point and increase h until this fixed-point becomes unstable. We 
then hold h fixed and allow the system to evolve to a new fixed-point characterized by a higher m*(h). Holding h 
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constant during the evolution of the system amounts to the assumption that the applied held varies infinitely slowly as 
compared with the internal relaxational processes of the system. The above process is repeated until all fixed-points 
in increasing applied field are determined. The trajectory of the magnetizations m*(h) of these fixed-points forms the 
lower half of the hysteresis loop. The upper half of the hysteresis loop is obtained similarly by starting with the fixed- 
point with m* = 1 and decreasing the field in smallest steps to obtain the sequence of fixed-points up to m* = — 1. 
If the size of the sample N is sufficiently large so that finite size effects can be neglected, the magnetization m* u (h) 
on the upper half of the hysteresis loop is related to the magnetization m*(h) on the lower half by the theoretical 
symmetry relation m*(/i) = —m*(—h). Hysteresis loops for the undiluted case (c=l) have been studied numerically 
on (i-dimensional regular lattices, and by an exact solution on a Bethe lattice H G3. In the following we extend the 
exact solution on the Bethe lattice to c < 1 and apply it to rock magnetism. 



III. HYSTERESIS ON A BETHE LATTICE 



A Bethe lattice is an infinite-size branching tree of coordination number z. It may be visualized as the deep interior 
of a large branching tree (Cayley tree) of the same coordination number. Figure (1) shows a small ( 4 levels) Cayley 
tree with z = 3 drawn such that the lattice points at the bottom row (level 0) are at the surface of the tree, and the 
point at the top (level 3) is at the root (center) of the tree. All lattice points except on the surface have z nearest 
neighbors. A lattice point on the surface has only one neighbor. Analysis of hysteresis on a Cayley tree is simpler 
because there are no closed loops on the lattice. However most of the lattice points lie on the surface and therefore 
special care has to be exercised to ensure that the results apply to the deep interior of the tree and are insensitive to 
conditions on the surface. We adopt two separate methods to eliminate the surface effects. In our analytic calculations 
we use recursion relations that take us from the surface towards the interior of the tree level by level. Fixed-points 
of these recursion relations arc insensitive to the random-fields on the surface. In numerical simulations we employ 
a different strategy. We perform simulations on a surfaceless random graph of coordination number z. The two 
methods of eliminating surface effects are equivalent and therefore our theoretical results match our simulation results 
perfectly. 

Consider a tree whose sites are randomly occupied by an Ising spin with probability c. If c is slightly less than 
unity the lattice acquires holes (regions without spins). For values of c closer to zero it breaks up into disjointed 
clusters of spins. It is not immediately obvious that the earlier method of recursion relations [i| starting from the 
surface of a compact lattice (c = 1) and moving towards its interior is still useful i.e. if it would yield fixed-points 
in spite of enhanced surface effects. But we find that the method works over the entire range < c < 1. We start 
with a diluted lattice and h — — oo so that all spins are down initially. The field is slowly ramped up from — co to h 
and we ask what fraction of spins are up at h. We choose a site at random in the deep interior of the tree and call 
it the central site. The probability that the central site is occupied by a spin is equal to c. We now calculate the 
probability that it is up as well. Each nearest neighbor of the central site if it is occupied by a spin forms the vertex 
of a subtree. The subtrees do not interact with each other except through the central site. Therefore the evolution on 
z subtrees meeting at the central site is independent of each other as long as the central site does not flip up from its 
initial state. It is an important property of the model that the order in which the spins flip up does not influence the 
fixed-point at h so we may assume that the central site is the last site to flip up at h if it flips up at all. Before 
we can say whether the central sites may flip up or not, we need to know the conditional probability cP*(h) that a 
nearest neighbor of the central site is up given that the central site is down at applied field h. 

P*(h) is the fixed point of the recursion equation, 



i=0 



- 1 

i 



(i-cf 



E 

J=0 



■ l-i 

j 



{cP'" 1 ^)} {c-cP l ~ l {h)} 



I z-l-i-j 



Pj(i;h) 



(5) 



The above equation is written on the assumption that all spins in the system were down before being exposed to 
a field h. Spins on the surface had the first chance to flip up at h, then spins on level 1 and so on. In this process of 
organized relaxation, spins up to level I — 1 have been relaxed and we are at the point of calculating the probability 
that a spin at level / (say at site x) flips up while spins at level I + 1 are still down. With this explanation, equation 
(5) is easily understood after various symbols are defined. P l {h) is the conditional probability that the spin at site x 
at level I is up given that its nearest neighbor at level I + 1 is down. Similarly, P i_1 (/i) is the conditional probability 
that a site at level I — 1 is up given that its nearest neighbor at level I (site x) is down. The site x has one neighbor 
at level I + 1 and z — 1 neighbors at level I — 1. These neighbors can be occupied by spins with probability c or 
unoccupied with probability 1 — c. The neighbors at level I — 1, if occupied, could be up or down independently of 
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each other; Pj(i; h) is the probability that site x has sufficient quenched field to flip up at h if j neighbors arc up, i 
neighbors are unoccupied by spins, and consequently z — j — i neighbors are down. 



Pj(i;h) 



</>(hi)dhi; 4>(hi) = 



1 



>(z-2j-i)J-h V27Rj2 

The probability that a randomly chosen site on the lattice (the central site) is occupied and up is, 



(6) 



P(h)=c^(. {l-c} 



J2( . ){cp*{h)y{c-cp*{h)y-^p^,h) 

3=0 V J / 



(7) 



The magnetization per spin on the lower and upper half of the major hysteresis loop is given by 



mUh) = 2p(h) 



m *u( h ) = -m*i{-h) 



(8) 



If we reverse the applied field before completing the lower half of the major loop, we generate a minor hysteresis 
loop. First reversal of the field generates the upper half of the minor loop, and a second reversal generates the lower 
half. When the field on second reversal reaches the point where the first reversal was made, the lower half of the 
minor loop meets the starting point of the upper half. In other words, the minor loop closes upon itself at the point 
where it started. This property of the random-field Ising model is known as the return point memory. Consider the 
upper half of the minor loop. Suppose the applied field is reversed from h to h d (h d < h). We need to calculate the 
probability that an occupied site, say site x, that is up at h turns down at h d . When site x turns up at h, the field on 
its nearest neighbors increases by an amount 2 J. This may cause some neighbors to turn up as well. Each neighbor 
that turns up increases the field on site x by an amount 2 J. Therefore on reversing the field, site x can turn down 
only after all neighbors which turned up after it have turned down. The probability D*{h d ) that an occupied nearest 
neighbor of site x that was down before site x turned up at h is again down at h d is determined by the fixed point of 
the following recursion relation. 



D*(h d ) = cJ2{\ 1 )( 



z — l~i 

E 

3=0 



z — 1 — i 



{cP-(h)}'{c-cP'(l,)r- l -'- l {l-p^i(i;h)} 



(*;>-)* 



E 

3=0 



z — 1 — i 
j 



j'.rr>*fA d u z - 1 -i- i 



{cP*(h)y{D*(h d )} 



{p j+1 (i;h) - p j+1 (i;h d )} 



(9) 



Given an occupied site x that is up at h, the first sum above gives the conditional probability that an occupied 
nearest neighbor of x remains down at h after site x has turned up. The second sum takes into account the situation 
that the nearest neighbor in question turns up at h after site x turns up but turns down at h d . 

The fraction of occupied sites that turn down at h d is given by 



l(h d )=cJ2 (* (l-c) 



i=0 



E 

3=0 



z — I 

j 



{cP*(h)y{D*(h d )}'- i ^ip j (i,h)- Pj (i,h' t )} 



(10) 



The magnetization on the upper return loop and in the range of the applied field h — 2J<h d <his given by 



m(h d ) = 2{p(h) - q(h d )} - c (11) 

At h d = h — 2J, all occupied neighbors of the central site that flipped up because the central site flipped up at 
h would flip down but the central site would stay up. If the applied field decreases further, the central site would 
turn down before any of its occupied nearest neighbors. This means that at h — 2J the system arrives at some point 
on the upper half of the major hysteresis loop, and moves on it upon further decrease in the applied field . The 
magnetization for h d < h — 2 J is given by 



6 



m{h d ) = 2p(h d ) - c 



(12) 



Where, 



E 

3=0 



{cP*(ft d )H{(i - cP*{h d ))y-^ P] {z- h d ) 



(13) 



and P*(h d ) is given by the fixed point of the following recursion relation. 



E 

3=0 



z — 1 — i 



{cP , - 1 (/i d )} J '{c-cP'- 1 (/i < ')} 



l ~ 1 lh d \\ z - 1 ^ i ^i 



p j+1 {i;h d ) 



(14) 



The lower half of the return loop is obtained by reversing the field from h d to h u (h u > h d ). If h d < h — 2 J, the 
lower half of the minor loop starts from the major loop, and therefore it is related by symmetry to the upper half 
of the return loop that has been obtained above. We only need to consider the case h d > h — 2 J. In this case, the 
magnetization on the lower half of the return loop may be written as, 

m(h u ) = 2{p(h) - q(h d ) + p'(h u )} - c (15) 
Where p'(h u ) is the probability that an occupied site that is up at h, down at h d , turns up again at h u . 



E 

3=0 



J 



{u*(h u )y{D*(h d )y-^{ P] ( t , h u ) - Pj (i, h d )} 



(16) 



Here U*{h u ) is the conditional probability that an occupied nearest neighbor of a site x turns up before site x turns 
up on the lower return loop. It is determined by the fixed point of the following equation. 



U*(h u ) = cP*(h)-cJ2 

i— 
z-l 

+^E 



i=0 



z-l 



i 

z-l 



(1-c)' 

(i- c y 



3=0 V 3 

E ( z ' 1 ' 1 ) 

3=0 ^ 3 



(17) 



The rationale for the equation (17) is as follows. Given an occupied site x that is down at h d , the first two terms 
on the right hand side account for the probability that an occupied nearest neighbor of site x is up at h u > h d . Note 
that the neighbor in question must have been up at h in order to be up at h d , and if it is already up at h d , then it will 
remain up on the entire lower half of the return loop, i.e. at h u > h d . The last term gives the probability that the 
nearest neighbor was down at h d , but turned up on the lower return loop before site x turned up. It can be verified 
that the lower return loop meets the lower major loop at h u = h and merges with it for h u > h, thus proving the 
property of return point memory. 

The method of calculating the minor loop may be extended to obtain a series of smaller minor loops nested within 
the minor loop obtained above. The key point is that whenever the applied field is reversed, a site x may flip only 
after all neighbors of site x which flipped in response to the flipping of site x on the immediately preceding sector 
have flipped back. The neighbors of site x which did not flip on the preceding sector in response to the flipping of 
site x do not flip in the reversed field before site x has flipped. We have obtained above expression for the return 
loop when the applied field is reversed from h on the lower major loop to h d (h — 2 J < h d < h), and reversed again 
from h d to h u {h u < h). When the applied field is reversed a third time from h u to h dd (h dd < h u ), expressions for 
the magnetization on the nested return loop follow the same structure as the one on the trajectory from h to h d . 



7 



The analytic results obtained above are depicted in Figure (2) for a = 1.7 and z = 4 for two concentrations of 
magnetic minerals; (i) c=l (no dilution), and (ii) c=0.8 . The fixed point of equation (5) is evaluated numerically for 
a sufficiently large and negative applied field h, and used in equations (7) and (8) to obtain the magnetization at the 
start of the lower half of the hysteresis loop (m* pa — c). The applied field h is then increased in small steps, and the 
process is repeated to determine to* (ft.) along the lower half of the hysteresis loop. Similarly, minor loops are obtained 
using equations (9) to (17). The upper half of the hysteresis loop is obtained by the relation m*(/i) = —m?(—h). 
Hysteresis on a z > 4 lattice is qualitatively different from the case of z — 2 and z = 3. For z > 4 there exists a 
non-equilibrium critical point (a = <r Cl h = h c ) on the lower half of the hysteresis loop, and another symmetrically 
placed critical point on the upper half. At these non-equilibrium critical points the response of the system to a slowly 
varying driving field is singular. Critical points do not exist on lattices with z = 2,3. Let us focus on the critical 
point on the lower half of the hysteresis loop. For z = 4 and c = 1, we have a c pa 1.78. For z > 4, a c increases 
with increasing z. For a < a c , the two halves of the hysteresis loop have first order jumps in the magnetization as 
shown in Figure (2). For c < 1, there is extra disorder in the system on account of the dilution of magnetic sites. 
This causes the first order jumps at a = 1.78 for c = 1 to vanish at c = 0.8, and the major hysteresis loop becomes 
smooth as seen in Figure (2). We have also shown two minor hysteresis loops for c = 0.8. If c < 1 but sufficiently 
large to form spanning clusters of occupied sites on the lattice, the qualitative behavior of hysteresis is similar to 
that on the undiluted lattice with c = 1 but with an effectively reduced value of a c . These results are verified by 
numerical simulations of the model for the corresponding choices of the parameters z, J, c, and a. As stated earlier, 
we performed numerical simulations on random graphs to eliminate surface effects. A random graph of N sites has 
no surface, but the price we pay is that it has some loops. However, for almost all sites in the graph, the local 
connectivity up to a distance of log^ z _i^N is similar to the one in the deep interior of the branching tree. Thus the 
the theoretical results depicted in Figure (2) to (4) perfectly matched the corresponding results from the simulations 
of the model on random graphs. Indeed, the theoretical and simulation results are indistinguishable on the scale of 
the figures. 



IV. APPLICATION TO ROCK MAGNETISM 



There are a variety of models in the literature that explain the shape of hysteresis loops in magnetic rock samples. 
Our objective is not merely to add to this list of models but also to simplify the situation. Our point is that extant 
models of ferromagnetic hysteresis also explain hysteresis in magnetic rocks if we add the ingredient of dilution to 
them. We have chosen the random field Ising model of hysteresis in ferromagncts at zero temperature and in the limit 
of zero frequency of the driving field. The advantage of this model is that it can be solved exactly on a Bethe lattice 
of an arbitrary coordination number z [ToT ]. The model has only a few parameters; the coordination number z of the 
lattice, the interaction energy J that aligns nearest neighbor spins parallel to each other, and a parameter a that 
tends to disorder the system. In spite of its simplicity, this model explains a number of experimental observations 
regarding return point memory, Barkhausen noise, non-equilibrium critical points on hysteresis loops, and universal 
behavior in the vicinity of these points [HQ. AH we have done is diluted this model, i.e. only a fraction c of the sites 
are occupied by magnetic entities. We have solved the dilute model exactly. Drawing upon the empirical observation 
that the magnetic minerals in a rock account for only a few percent of its mass, one may ask if the predicted hysteretic 
behavior of our model in the range 0.01 < c < 0.1 is reasonably close to the observed behavior of magnetic rocks. It 
indeed appears to be the case considering the minimal number of free parameters in our model. 

Figure (3) shows three hysteresis loops for c = 0.1; and a = 0.05,0.2, and 0.5 respectively. We have set J = 1 for 
convenience. Let us focus on a = 0.05 first. As the applied field h increases from h = — oo the magnetization stays 
at its saturation value m = —0.1 until h pa — 3a. It rises sharply from to sa —0.1 to m pa 0.03 around h sa 0, and from 
to pa 0.03 to to pa 0.1 around h sa 1. The rise around h = is due to isolated spins in the system. The quenched 
field has a Gaussian distribution with average value zero and standard deviation a. Therefore some of the isolated 
sites turn up at h sa — 3a, and almost all of them are up at h pa 3a. The fraction of isolated sites is c(l — c) 4 . This is 
approximately equal to 0.034. Occupied sites in clusters of size larger than unity account for the remaining fraction 
0.066. After all isolated sites have turned up the magnetization per site increases to —0.1 + 2 x 0.34 = .032. It stays 
at this plateau value until the applied field increases further by an amount 2 J — 3a enabling some surface sites of a 
connected cluster of down spins (dangling bonds) to turn up. Nearly all sites turn up at h = 2 J + 3a and we get the 
saturated magnetization to = .1 as shown in figure (3). The same mechanism applies to hysteresis loops for other 
values of a. With increasing a the sharp corners tend to get more rounded and the middle (pinched) portion of the loop 
broadens. We get a gentle pinched loop for a = 0.2. For a = 0.5 the constriction in the middle disappears completely 
and we get a normal ferromagnetic type of hysteresis loop as shown in figure (3). Thus our model reproduces various 
shapes of hysteresis loops observed in magnetic rocks with a minimum number of parameters and without invoking 
different mechanisms for different shapes. Figure 4 (see caption) shows another way of representing the same data 
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as shown in figure 3. In this representation, the constriction in the middle of the hysteresis loop appears as a more 
pronounced dip in the middle of the plot and therefore it may be of some value in the analysis of experimental data. 

Can we deduce the magnetic composition of the rock from the shape of its hysteresis loop? This is not possible 
in general because the hysteresis loop contains the information in an integrated form. Information is irreversibly 
lost in the process of integration and can not be retrieved. There is no unique way to obtain the components of a 
sum from the sum! However, our model provides a few guiding principles. Ferromagnetic shapes of hysteresis loops 
indicate a larger random-field disorder a in the system. On the other hand a wasp-waisted loop especially with a 
long and narrow waist and sharp bends indicates relatively small disorder a. In this case the connectivity of spins 
at the surface of clusters has the larger effect on the shape of hysteresis loops. It should be possible to deduce the 
magnitude of interaction J from the locations of sharp turns in the magnetization along the applied field. Similarly, 
it should be possible to deduce the relative size of clusters from the plateaus in the magnetization. The distribution 
of cluster sizes and the plateaus of the magnetizations are known exactly in our model in terms of the two parameters 
(J and a) of the model. However, more complex rock samples may not be adequately characterized by a simple 
two-parameter Hamiltonian. Also, there may be systems whose frozen disorder is not represented adequately b 
a Gaussian distribution. We have focused on the Gaussian distribution because it applies to many systems |2|, 
and its use is justified by the central limit theorem. The effect of the shape of random-field distribution on critical 
hysteresis in the random- field Ising model (c = 1) has been examined by Liu and Dahmen (28|. They find that the 
Lorentzian and parabolic distributions of random fields yield the same critical exponents in three dimensions as the 
Gaussian random fields. This is not entirely surprising because a renormalization group analysis in 6 — e dimensions 
[29j shows that for random-field distributions with a single maximum, it is only the curvature of the distribution at 
its maximum that determines the critical exponents. Therefore it is reasonable to use the Gaussian distribution as 
a common representative of smooth single-peaked distributions which all give the same critical exponents that agree 
with experiments. In the same vein, we have used the Gaussian distribution to understand the shapes of hysteresis 
loops in the strongly diluted limit of the random-field model. Before ending, we also wish to mention that wasp- 
waisted loops are not a property of magnetic rock samples only. Such loops are also seen in random magnets with 
the order parameter having a continuous symmetry shape memory alloys l30l |3ll . and martensites [33 ] . Indeed 



the physics behind hysteresis in the random field Blume-Emery-Griffiths model [32| is very similar to that discussed 
here in the case of the dilute random-field Ising model. 
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FIG. 1: (Color online) A Cayley tree with with 4 levels (I = 0, 1, 2, 3) and coordination number z — 3. Each node is connected 
to z nearest neighbors except the nodes at the surface (Z = 0) which have only one neighbor. The deep interior of the tree 
where surface effects can be neglected is known as the Bethe lattice. 



-1.5 -0.75 0.75 1.5 

h 

FIG. 2: (Color online) Hysteresis loops on a z — 4 lattice for a = 1.7: The larger loop (red/dark) with saturation at m — ±1 
and first-order jumps in the magnetization is for c = 1 (undiluted lattice); the smaller loop (green/grey) with saturation at 
m = ±0.8 is for c = 0.8; Two minor loops within the smaller loop are also shown which are obtained by making excursions 
from h — 0.9 to h = —1.1 and h — —0.85 respectively. 
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FIG. 4: (Color online) This is a different representation of figure (3). For a fixed magnetization m along the ic-axis, the y-axis 
shows the difference of corresponding applied fields on the upper and lower halves of the hysteresis loop. The red line (nearly 
horizontal thin line) corresponds to a = .5, the blue line (thick line with two nearly vertical jumps) corresponds to a — 0.05. 
The green line (intermediate dotted line) corresponds to a — 0.2. 



